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Abstract: We present a new perturbative approach to the study of 

signal-noise interactions in nonlinear optical fibers. The approach is based 
on the hydrodynamic formulation of the nonlinear Schrddinger equation 
that governs the propagation of light in the fiber. Our method is discussed in 
general and is developed in more details for some special cases, namely the 
small-dispersion regime, the continuous-wave (CW) signal and the solitonic 
pulse. The accuracy of the approach is numerically tested in the CW case. 
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1. Introduction 

In fiber-optic communication systems, optical amplifiers are usually employed to periodically 
compensate for fiber loss. As a side effect of amplification, amplified spontaneous emission 
(ASE) noise is added to the optical signal at each amplification stage. Both signal and noise 
propagate through the optical fiber according to the nonlinear Schrodinger equation (NLSE), 
which accounts for attenuation, dispersion, and Kerr nonlinearity. In the linear regime (i.e., at 
low power), the ASE noise is not affected by propagation through the fiber and can be modeled 
as additive white Gaussian noise (AWGN) at the end of the optical link. On the other hand, 
when nonlinear effects are not negligible, signal and noise interact during propagation. This 
interaction manifests, for instance, in the parametric gain and modulation instability effects H] 
and in the emergence of nonlinear phase noise |l2l. As a result, noise becomes, in general, col¬ 
ored, non-stationary, and non-Gaussian. The knowledge of the statistics of the noisy signal at 
the output of the link is extremely important to design efficient optical systems, evaluate their 
performance, and establish the ultimate channel capacity. Eor instance, optimum detection (in 
a maximum likelihood or maximum a-posteriori probability sense) is based on the knowledge 
of the conditional distribution of the output signal given the input Q. Unfortunately, this dis¬ 
tribution is, in general, unknown. Exact statistical models in the absence of fiber dispersion 
have been studied in lldUbl. In the presence of dispersion, approximate models based on per¬ 
turbation methods have been investigated in Some alternative (non-regular) perturbation 

methods have been also proposed to describe the noise as a non-additive perturbation of the 
noise-free solution. Eor instance, the combined regular-logarithmic perturbation (CRLP) ifTOl 
provides an accurate description of the non-Gaussian distribution through a non-linear combi¬ 
nation of Gaussian variables, accounting for both parametric gain and nonlinear phase noise. 
Moreover, in |[3l it is shown that a polar-Gaussian metric may provide a better performance 
than a Cartesian-Gaussian metric for maximum likelihood sequence detection. This suggests 
that a perturbation approach applied to a polar representation (e.g., amplitude and phase) of the 
signal may be a good alternative to obtain a simple but effective statistical model. 

A similar problem arises in the different context of high-power fiber lasers, in which an 
incoherent or partially-coherent CW light is subject to spectral broadening during propagation 
through an optical fiber EHil. Also in this context, the important role of the interaction 
between fiber dispersion and nonlinearity and the suitability of a polar representation of the 
signal have been widely recognized. 



In this paper we develop a novel approach to the problems outlined above, which is based 
on the hydrodynamic formulation of the NLSE. Such formulation, originally introduced by E. 
Madelung in a quantum mechanical framework ifTsll and subsequently applied also to optical 
fibers (e.g, in ifTMTgl ). consists of a system of compressible, isothermal, Euler-like equations 
with a typical “quantum pressure” term, also known as Bohm potential. In such system the 
hydrodynamic variables have the physical meaning of power and angular frequency (playing 
the role of the fluid density and velocity, respectively). When the dispersion is small compared 
to the nonlinearity, the quantum pressure term can be formally neglected and we are left with 
a much simpler hydrodynamic system which is similar in form to the so called “shallow-water 
equations”. Such a semiclassical approximation has been employed, for instance, in the study 
of self-similar propagation of high-power pulses in optical fiber amplifiers Il20l - l23l . It must be 
noticed, however, that in this way the dynamics loses its dispersive character, which may lead 
to formation of shocks or even to ill-posedness of the mathematical problem ll24ll . Both the full 
and the approximated hydrodynamic equations are then used to study the propagation of an 
input signal affected by the ASE noise. Assuming the mean power 2 ( 7 ^ of the injected noise to 
be small with respect to the typical power of the signal, we perform a perturbation expansion at 
first-order in the small parameter a. The leading order describes the fully deterministic propa¬ 
gation of the signal while the first-order obeys a linear system with stochastic input. Assuming 
the input ASE to be described by a band-limited white gaussian process, it turns out that, at 
first order in a, the input power and frequency are still described by a joint gaussian process. 
Hence, the linearity of the stochastic system implies that the statistics remain gaussian after the 
propagation, which simplifies the computation of the output statistics. 

The paper is organized as follows. The NLSE is recast in hydrodynamic form in Section |2l 
where also the small-dispersion regime is considered. In Sec. |3l the statistics of the hydrody¬ 
namic variables are computed up to terms of order <7^. The perturbative expansion is introduced 
in general in Sec. |4]and is applied to the particular regimes in the subsequent sections |5](small- 
dispersion regime), | 6 ](CW signal) and|2](solitonic pulse). In Sec. [ 8 ]the perturbative method is 
used to get equations that govern the evolution of the expected value of the solution. Einally, 
Sec.|9]is devoted to numerical experiments, in the CW case, that are aimed to test the accuracy 
of the method. 


2. Hydrodynamic formulation of the NLSE 


We assume, as usual, that the electromagnetic field propagates along the fiber according to the 
NLSE ED 


= 


, I ,2 .« 

-yMu + rlwl U-I—U, 


( 1 ) 


where the subscripts denote derivation with respect to z and t. Here, z > 0 is a longitudinal 
coordinate of the fiber axis, f e R is the time coordinate corresponding to a time frame moving 
with the signal group velocity, j 32 G R is the group velocity dispersion parameter and 7 > 0 
and a > 0 are, respectively, the nonlinear refraction and attenuation parameters. The variable 
u = u{z,t) is the envelope of the modulation of the electromagnetic field along the fiber. The 
dispersion parameter j 32 can be either positive (normal dispersion), or negative (anomalous 
dispersion), giving the NLSE, Eq. O, a defocusing or focusing character, respectively. 

We apply to Eq. ([T]i the transformation 


M = V^e''^, v = 0r, 


( 2 ) 


where the real quantities p > 0, (j) and v', which are functions of z and t, are, respectively, 
power, phase, and angular frequency of the signal u. Substituting Eq. ^ into Eq. ([TJ, and 


writing separately the equations for the real and imaginary part, we obtain the “hydrodynamic” 
formulation of the NLSE: 


' Pz + P2ipv)t = -ap, 

Apart from the exchange of space and time coordinates, the presence of the dissipative term 
— ap and, of course, different physical constants, system Q is the hydrodynamic form of 
Schrodinger equation (for the nonlinear potential V = JP) that was introduced by E. Madelung 
in 1926 Ha. They have the form of an Euler system of conservation laws for an isothermal, 
compressible gas with a “quantum pressure” term P 2 Q, where 

^ {y/p)tt ip,f-2ppt, 

^ 2VP 8p2 • ^ ^ 

By considering a reference time Tq and power Pq (determined by the duration and mean power 
of a pulse injected in the fiber), it is known that the dimensionless parameter 


£ 


1 [M 

To Y yPo ’ 


(5) 


measures the relative strength of the dispersion term ll^ (it plays the role of a scaled Planck 
constant for the NLSE ([Til)- By using dimensionless variables, it is possible to show that the 
error in neglecting the term in Eq. (O is of order e^. Then, for £ <C 1 we can formally 
write the approximated system 


f Pz + p 2 ipv), = -ap, 

{ 7 \ (6) 

[>■-.+{y'- +7p),=«- 

Eor rigorous results about the limit £ —>^ 0 of system (|3) in the defocusing case, one can refer 
to El. 


Remark 2.1 Condition £ <C 1 can be interpreted as a condition on the typical frequency 




(7) 


which is generally not satisfied in typical optical fiber communication systems.Thus, while the 
approximated system (jh]) can be used for specific applications (e.g., high-power fiber lasers or 
optical amplifiers), we will employ the complete system (H to obtain some results of more 
general validity, suitable also for optical communications. 

Remark 2.2 Equation (|6]i with a = 0 and j32 > 0 is a hyperbolic system of conservation laws 
(indeed, it is formally equivalent to the well known equations describing shallow waters). Eor 
j 32 < 0 such system is not hyperbolic, which implies that the input-value problem is not well 
posed in general and leads to unphysical results E6l . This can be clearly seen in the simulations 
reported in Sec.|9l where the small-dispersion approximation applied to a fiber with anomalous 
dispersion predicts a non-physical, exponential growth. 


System (|3), or its small-dispersion approximation, Eq. (|6]), will be supplemented with input 
conditions (i.e. conditions at z = 0) of stochastic nature, which will be described in next section. 








3. Stochastic input data 

We assume that an optical amplifier injects in the fiber, at z = 0, an input signal of the form 


M(0,f) = Mo(0 + O’Z(f) 


( 8 ) 


where uo{t) is the deterministic part of the input signal, and aZ{t) is the stochastic part, that 
is the ASE noise. We assume that Z(f) = Zi(f) + iZ 2 {t) is a gaussian, band-limited, complex 
white noise lIZTll . i.e. a stochastic process such that: 

i) for every fixed f G M, Zi (f) and Z 2 (f) have normal distribution with mean 0 and variance 
1 (so that the variance of the noise is 2(7^), and are independent; 

ii) the autocorrelation function is 



(9) 


(where (•) denotes the expectation), corresponding to a power spectral density 



( 10 ) 


where Wc > 0 is the cutoff frequency. 

Note that 2(7^ is the mean power of the noise and in the following we will assume 

( 7 ^ < Po. 


( 11 ) 


3.1. Statistics of the input power and phase 

The input condition, Eq. (|8ll, has now to be translated into an input condition for the hydrody¬ 
namic system, Eq. I©. This will be done at first order in <7 (i.e. up to terms of order that are 
small compared to the typical pulse power Pq). 

As far as the power is concerned, we can write 


\uq -f- cZ] — [mqI -f 20”Re(MQZ) -f ^(c7^). 


Hence, up to 

|M(0,f)P = p(0,f) ft! Po(f) + C7R(f), 

where, according to Eq. 0, the deterministic part of the input has been decomposed as 

“0 = VPo^“^°- 


( 12 ) 


(13) 


The stochastic perturbation of the input power. 


R{t)=2Re{uo{t)Z{t)), 

is therefore a gaussian process with 0 mean and (by Eq. (|9]l) autocorrelation 


(14) 


{R{t)R{s)) = Ay/po{t)Po{s) cos(^o(0 - (I>q{s)) sinc(®t.(f - s)). 


( 15 ) 


As far as the phase is concerned, we have 







( 16 ) 


Since (j) = arg(M), we have that, up to 

where the perturbation 


4>(f) = - 


\m[uQ{t)Z{t)) 

Poit) 


is a gaussian process with 0 mean and autocorrelation 


^\ Re(Mo(0«o('S)) ■ , u \\ 


(17) 


(18) 


3.2. Statistics of the input frequency 

We can now compute the statistics of the input frequency 


^(0,?) = —0(o,f) Ri vo(f) + (jy(f), 


(19) 


where, vo(f) = 0o(f) and V{t) = By standard results on differentiated processes ll27l we 
have that V{t) is a gaussian process with 0 mean and autocorrelation function 


{V{t)V{s)) = ^{^{t)^{s)). 


dtds 


Then, by using Eq. CHl, we obtain 


(y(f)y(i)) = |®c Re — Re sinc'(®t(f — i)) 


Re 


sinc(®c(f — s)) — Re sinc"(®c(f — s))|, (20) 


where 




Po ypo 


and tJ' = 


VPo 


0 V 


ivo- 


A. 

2po 


We finally compute the cross-correlation between R and V. First of all, from Eqs. (01, (dll 
and (fTTl i. we immediately obtain 

(f?(f)4>(s)) = 2Im sinc(cUc(f — s)). (21) 

Then, using {R{t)V{s)) = ^(/?(f)<t>(s)), we have 
(R(f)y(s)) = 2Im sinc(cUc(f — s)) — 2a)tlm ^Mo(f)r?(s)^ sinc'(a)c(f — s)). (22) 


4. Perturbative treatment of the stochastic problem 

Let us rewrite here the stochastic input-value problem given by Eqs. (0, (fT2l i and ( [10 1: 


Pz + P2{pv)t = -ap, 


t'z + 



(23a) 


p{0,t) = Po{t) + aR{t) 
v(0,f) = vo(f) + (7y(f) 


(23b) 

















This problem will be treated perturbatively with respect to the small parameter a; then we 
expand the unknowns as 

p + V = + ^(c7^), 

where we assume p^^^ ~ ~ <7. Substituting into Eq. ( l2Tl i yields, at leading order, 




(24a) 


with input conditions 


'p(“)( 0 ,f) = po(f) 

v^°'>{ 0 ,t) = vo(t) 


(24b) 


(which, of course, means that at leading order we are left with the deterministic problem) and, 
at first order. 


.(1) 


+ + =-ap(') 


+ 


^ypW 


= 0 , 


' p^^'>iO,t) = aRit) 
v(^'(o,f) = c7y(f) 


where is the first-order approximation of Q (see Eq. (IHi), given by 

,(l) l/n(')N ^ 

4p(o) 


n\^l 1 nil) 1 /nl‘l\ 


(25a) 


(25b) 


(26) 


(linear in p (^^). Note that Eq. (l24li is a nonlinear deterministic problem (equivalent to the NLSE) 
while, once Eq. (l24l) has been solved, Eq. (|25] | is a linear stochastic problem. 

Remark 4.1 The stochastic problem introduced in this section refers to the case of signal and 
ASE noise injected at the input of a single span of fiber. The analysis can be easily extended to 
a multi-span optical system by following the same approach employed in Eol. In this case, at 
each amplification stage, an independent perturbation of the order of a (due to the ASE noise 
generated by the inline amplifier), with statistics as in Eqs. ® , (|20l) and is added to the 
(already noisy) Madelung variables p and v, which are then propagated through the following 
span of fiber according to Eq. (I23ab . 


5. The small-dispersion approximation 

Let us consider the small-dispersion approximation {Q k. 0) of the input-value problems (l24t 
and (|25]) . The leading order equation, Eq. (I24ab . becomes 


pP + P2{p^^'’v^°'’)t = -ap^^\ 

.•f + (|M"V + rp<"'),=o, 


( 27 ) 







that is the “shallow water” system (|6l). As already mentioned in Remark 12721 according to the 
standard theory of first-order systems II26II28I . we have that the system is hyperbolic if and only 
if & > 0. In this case the characteristic velocities are 

c±=j32v(‘’)±y&W, (28) 

and the characteristic curves are defined by the equations 

| = c±=j32v(‘')±^j327p(«). (29) 

Moreover, we obtain the equations 


d 

dz 




= -a 



on 



(30) 


for the Riemann variables 2y/p ^± For a = 0, the Riemann variables are constants 

along the characteristic lines, in which case they are Riemann invariants for the problem. The 
condition a = 0 can always be assumed on a short fiber stretch, provided that the nonlinear 
coefficient 7 is substituted with an effective coefficient 7 eff 

Remark 5.1 Although the characteristic equations, Eq. (l30b . do not lead to the analytic solu¬ 
tion of Eq. (I27I 1 (because they hold separately on the two characteristic curves, which depend on 
the solution itself), nevertheless the form Eq. (l30l l allows a very efficient numerical treatment 
of the problem (see e.g. ED). 

Assume now that is the solution to problem (IZTT i. and consider the small-dispersion 

approximation of the first-order equation (125 ab : 

This is a linear, non-autonomous, problem for the unknowns (p(^(, Eor ^2 > 0, its char¬ 
acteristic curves are fixed, and are the same as those of problem (IZTl i. i.e. they are still given by 
Eq. ( |29] |. The characteristic equations read as follows; 


^r±+j32v(‘’V± + 
dz 







dt 

on — =c±, 
dz 


(32) 


where r± = ± ^ y(i) are the Riemann variables for the first-order problem. 

6. The continuous-wave case 

The simplest case in which problems (l24b and (l25l l can be explicitly solved is that of constant 
deterministic input data po and vq. We shall also assume a = 0. Under these conditions, problem 
(l24b admits the constant solution 


p(°)(z,f) = po, 


v(°^(z,t) = i'o- 


(33) 










Such solution is physically meaningful and corresponds to a continuous-wave (CW) solution 
of the NLSE O (with a = 0) HO) 


M(2,f)=e-'(fcT+rPo).enwy^ 

(up to a constant phase). The first-order equation, Eq. (I25al) . becomes 

pi ^ + PivqPi ^ ^ = 0, 

^2 ^(1)) 


-f j32Vovi‘' + (yP^^'’ - ^P«= 0 


(34) 


(35) 


At this point we can notice that we can get rid of the terms containing vo by means of the change 
of variable t^t — p 2 Voz (indeed, vq just corresponds to a global frequency shift), which leads 
to 

pi'^+fepov'^' =0 




= 0 . 


(36) 


In order to make a comparison with the CRLP method ifTOl . we return to the the phase variable 
= V for which we easily obtain 


I pi‘’+j32Po0i/^ =0, 

Ap//) =0, 


with the stochastic input conditions 




(37a) 


(37b) 


Note that an arbitrary function of z was set to 0 in Eq. ( I37al ); this can be justified by looking at 
the Madelung equation for the phase. 


0z + y(0r)^ + 7P +(322 = 0, 

of whom Eq. ( I37al ) is the linearization. The comparison of Eqs. ( iTTl i with the CRLP model is not 
straightforward since, in the latter, a three-variable perturbation expansion is used to represents 
the optical signal. However, as shown in Qol, the CRLP expansion can be reduced to a simpler 
amplitude and phase representation by setting the quadrature component equal to zero. This 
is a reasonable approximation when nonlinearity dominates over dispersion. In this case, the 
equation governing the propagation of the CRLP amplitude and phase variables is analogous to 
Eq. (I37al) . but for obvious differences (a factor of two dividing and multiplying the right-hand 
side of the first and second equation, respectively) due to the fact that Eq. (I37ab refers to a 
power and phase expansion. The Eourier transform of Eq. (I37al) is 


and the associated transfer matrix T (z, co), such that 


p(i) 

0(1) 


(z,®) =T{z,(o) 


aR\ 

a^J 




(39) 


can be explicitly computed and reads as follows; 

/ 

Tiz,co) = 


iiz'/k) 




P^ai^+^rPo 


sin(zv^) 






s(zv^) 


where 


k = k{co) = 


P2(0^ 

4 

.2 ^ 4ypo 


7Po 


(40) 


(41) 


This expression holds for jSz > 0, or jSz < 0 and co^ > The expression of T (z, (o) for jSz < 0 

and aP' < is easily obtained by using cos(«) = cosh(x) and sin(/x) = !sinh(x) (xreal). 

In the small-dispersion approximation, Eq. ( I37al ) drastically simplifies into 

I =0, 

1+ rp^^' = 0, 

and the corresponding transfer matrix, for jSz > 0, becomes 


(42) 


T{z.co) = 


cos(z®A/j32Po7) y^^®sin(z®A/i32Po7)\ 

_Z_ ^ cin -7/'ll. r*rvo ( -r/'il. / nA / 


^ sin(z®A/^2Po7) cos(z®a/^2Po7) 

while the expression for j32 < 0 can be obtained as explained above. 




(43) 


6.1. PSD matrix for the continuous wave 

The transfer matrix T{z,0)) allows to obtain the power spectral density (PSD) matrix 


G{z,co) = J 




\{<^^^\zf)p^^\zf + x)) (0(')(z,f)0^^H2T + T)),' ^ 

from its input value G(0, ®), through 

G(z,®) = 7’(z,®)G(0,®)r*(z,®), 




(44) 


(45) 


where Tf (co) denotes the adjoint matrix (which is just the transpose, in this case). The input 
PSD matrix 


G(0,®) = j 


' f{R(t)R{t + T)) (P(f)4>(f + T)) 

U^{t)R{t + T)) (<I>(f)0(f + T)) 


dr. 


(46) 


can be easily obtained from Eqs. (fTST i, ( fTSl l and (|2T]i, with mq = v^’ 0o = 0 and vq = 0, and 
reads as follows: 

na^ ( 4po 0 \ f (o 


G(0,(o) = 


[0 ^ 


1 rect 


V2®c 


( 47 ) 


Of course, in the “white noise limit” cOc —>-o° the rect function tends to 1. 















7. Perturbation of a soliton 


Let us now consider the one-parameter family of solitonic solutions of the focusing NLSE (i.e. 
Eq. ([T]i, with p 2 <0 and a = 0): 


4^)0 = d ■^t7sech(T7f)e'^'’^^/^, 


(48) 


where 77 > 0 is a parameter that determines the soliton width. A deterministic signal of the form 
( |48] | is a z-independent solution of Eq. ( I24ab with 


p{ 0 ) =772 sech^(77f), = 

7 2 




= 0 . 


(49) 


Substituting p and given by Eq. ( |49] | in Eq. (125al l we obtain an equation for the pertur¬ 
bation which, in the power-phase representation, takes the reads as follows; 


p_(i) _ Pi'H fsech^{r]t)d^dt = 0, 

'Y ^ / / 


0 : 


_|_ —ypW - 


( 1 ) 


(50) 


JL(_Pl_] 

4772 Vsech^fTif) A 


= 0 . 


2"' 4772 Vsech2(77f)/ 

Analogous results are obtained for the dark soli tons in the defocusing case (P 2 > 0): 

u{z,t) = W ^ 77tanh(77 7)e^'^’^^^, 


corresponding to 


= ^ri^tanh^{rit), = n5{t) - P 2 'n^z, v^^^ = n5'{t). 


(51) 


(52) 


Substituting and given by Eq. (0^ in Eq. (I25al) . and disregarding irrelevant terms 
coming from the delta function, we obtain the following equation for the perturbation, that we 
write directly in the power-phase representation; 


pi'^ -f ^^(tanh2(77f)^/'^)^ = 0 , 




csch 2 ( 77 f)sech 2 ( 77 f) \ ( 1 ^ 7 / p/^^ n _ 

FP 477Htanh2(77 7)A 


(53) 


2tanh2(77 7) 


It may be interesting to look at an approximation up to terms of order ae^. This amounts to 
keeping the leading-order problem, Eqs. dull, unchanged and neglecting the term containing 
(which is of order ae^) in the perturbed problem, Eqs. (l25T l. In the solitonic case, this 
means that the leading-order solution are still given by Eqs. dmi and (|48] |. but the second 
equation in both systems (l50l l and (l5^ is replaced by -f 7 p(^l = 0. Hence, in the focusing 
case we obtain the second order equation 

Pij^+i3|j?^(sech2(77f)p/'^), =0, 


( 54 ) 













which is of elliptic type. In the defocusing case, instead, we obtain the wave equation 


-i32J7^(tanh2(’?0p/‘’); = 0- (55) 

The elliptic character of Eq. (l54l i is another manifestation of the bad behaviour of the small dis¬ 
persion approximation in the focusing case, since it corresponds to the non-hyperbolic character 
of the first-order system for or (p('),v'(*)) (and see Remark |2^ . 


8. Equations for the expected value 

Let us finally consider the problem of obtaining equations for the expected values (p) and 
(v'). Since (R) = (V) = 0, taking the expectation of both sides of the first-order problem, Eq. 
(l25l l. yields a linear system in (p^'^) and (v(*)) with vanishing input data, which has the unique 
solution (p(^^) = (v(^^) = 0. 

Then, in order to see the effects of the noise on (p) and (v), we have to look at the next order 
in a: 

(p) = p^°^ -f (p(^^) -I- ^(cr^), (v) = -I- (v(^() -b 

In fact, by expanding the input Madelung variables p(0,f) and v(0,f), we see that the second- 
order system has a non-trivial input; 

(p(2))(0,f)=2cj2, (v(2))(0,f)=0, (56) 

(note that p(^)(0,f)/c7^ = |Z(f)p has a (2) distribution). The second-order equation for (p,v) 
is easily obtained from Eq. ( I23ab and reads as follows: 


pi^' + = -ap 


( 2 ) 


„( 2 ) 




= 0 , 


where is the second-order approximation of Q, given by 

_ 2p(‘’)p(2) - (p{l))2 /pW X l/p/^^x^ ^ 1 ^(p('i)^' 


0 ( 2 ) = 


(57) 
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Taking the expectation of Eq. (l57l) . and recalling that p(*’) is a deterministic quantity, we obtain 


(p(^))^-bj32(v(‘’\pp()-b(p(^'v(‘))-bp(‘’)(v(^^))^ = -a(pP'), 
(v(2))^ + [i32vW(vP)) + ^(v(')v(')) + 7{p^^^) + li2{Q^^% = 0, 


(59) 


which has to be supplemented with the input conditions ( |56] |. Note that Eq. ( l59b is a linear 
equation for the unknown (p(2)) which depends on the covariances between p(') and v('). 

This is an intricate problem in general, but it simplifies in particular situations. Eor example, 
let us consider the CW case, discussed in Sec.|6] In this case p^^'^ = po and = vq are two 
constants (and we can assume vq = 0 without loss of generality). Moreover, p(^) and v(*) are 
given by Eq. dJTl) and we have (p(^)p(*))^ = (p(*)v(^)), = = 0, since (p(*), v(*)) is a 

stationary process in the CW case. Then, from Eqs. ( l59b and ( fSST i, we obtain (in the power-phase 
representation) 

(p(2>)^H-j32po(0(^^)rt=O, 

R (60) 

(^{ 2 ))^ + ^(p{ 2 ))-A(p( 2 ))^^ = 0 . 








Hence, the equation for and in the CW case has the same form as Eq. dJTab . 

Supplemented with the input conditions (l5^ . it has the f-independent solution 

(pp))(z,f)=2(j2, (0(2))(o,r) = _27(j2^ (61) 

(up to an arbitrary phase), which simply represents the additional power due to the presence of 
noise and the corresponding nonlinear phase rotation due to the Kerr effect. 


9. Numerical results 

The accuracy of the proposed approach is tested in the CW case, for which exact analytical so¬ 
lutions for the PSD and pdf of the received signal have been derived. We consider the propaga¬ 
tion of a CW signal affected by band-limited AWGN noise through 50 km of lossless fiber. As a 
first test, we consider a non-zero dispersion-shifted fiber with GVD parameter P 2 = 5.1 ps^/km 
(normal dispersion) and nonlinear coefficient 7 = 1.27W^*km^^ a signal power of 10 mW, 
and a noise bandwidth of 50 GHz. The PSD matrix of the output power and phase variables is 
evaluated analytically from Eq. (l45l l. by using the transfer matrix either from the full first-order 
model, Eq. (l40l) . or from the small-dispersion approximation, Eq. (|4^ . Eor comparison, the 
same PSD matrix is evaluated by means of numerical simulations by employing the split-step 
Eourier method ll25l to propagate the noisy signal through the fiber, and the averaged peri- 
odogram to estimate the PSD. Eiguresllla) and[TJb) show the results obtained for an SNR of 20 
and 10 dB, respectively. In both cases, the significant deviation of the output PSD matrix from 
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Fig. 1. Normalized PSD matrix at the output of the non-zero dispersion fiber (normal dis¬ 
persion fh = 5.1ps^/km) for an SNR of (a) 20 dB and (b) 10 dB. 


the input one, Eq. (l47l i (which is diagonal and frequency-independent within the noise band¬ 
width), testifies the presence of a significant signal-noise interaction. At 20 dB of SNR, there is 
a perfect agreement between the full model and numerical simulations, while the approximated 
model, as expected, is accurate only at low frequency. On the other hand, at 10 dB of SNR, a 
small discrepancy between numerical simulations and the full model can be observed. This is 
due to the neglected higher-order terms in a which, at 10 dB of SNR, become relevant. 

In order to investigate the impact of fiber dispersion on the accuracy of the model, the same 
analysis is performed also for a standard fiber with a GVD parameter P 2 = —21.7ps2/km 
(anomalous dispersion), leaving all the other parameters unchanged. Eigures|2ja) and|2b) show 
the corresponding results for an SNR of 20 and 10 dB, respectively. 

In this case, the anomalous dispersion of the standard fiber is responsible for a significant 
amplification of the PSD of the power variable at some frequencies due to the presence of 
hyperbolic functions in Eq. ( l40b . This corresponds to an overall increase of the noise power due 
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Fig. 2. Normalized PSD matrix at the output of the standard fiber (anomalous dispersion 
j 82 = -21.7ps2/km) for an SNR of (a) 20 dB and (b) 10 dB. 


to signal-noise interaction. Moreover, the higher value of the modulus of the GVD parameter 
confines signal-noise interaction to a narrower bandwidth. Also in this case, the accuracy of the 
full model is confirmed, though the discrepancy at 10 dB of SNR is slightly more relevant than 
in the previous case because of the aforementioned amplification of the noise power. However, 
the approximated model is very inaccurate because of the higher dispersion. 

As a final test, we investigate the ability of the model to reproduce the non-Gaussian dis¬ 
tribution of the output signal. We consider the same CW power and the same kind of fibers 
considered in Figs. [Hand 121 an SNR of 20 dB, and reduce the noise bandwidth to 10 GHz 
(trying to keep all the noise within a bandwidth in which it is mostly affected by signal-noise 
interaction, as shown by the PSDs in the previous figures). Figures|3a) and|3jb) show the con¬ 
tour plots of the joint pdf of the real and imaginary parts of the (normalized and derotated) 
optical signal /^/po at the output of the fiber at an arbitrary time t (as, for a CW 

input signal, the process is stationary) for the non-zero dispersion-shifted fiber and the standard 
fiber, respectively. 




Fig. 3. Pdf of the optical signal at the output of the (a) non-zero dispersion-shifted fiber and 
(b) standard fiber. 

The solid curves are obtained according to the full first-order model, evaluating the PSD 
matrix of the power and phase variable as discussed before. Given the input conditions and 
the linearity of the propagation equations, the output variables are jointly Gaussian and their 
covariance matrix is evaluated by integrating the PSD matrix. The joint distribution of the real 
and imaginary parts is finally obtained through a change of variables. The dotted curves are oh- 

































tained by using the split-step Fourier method for the propagation of the noisy signal through the 
fiber, and applying the multicanonical Monte Carlo algorithm (see ll29ll and references therein) 
to estimate the pdf with uniform accuracy down to low values. In both cases, there is a good 
agreement between theory and simulations in the modal region of the pdfs, while some dis¬ 
crepancy appears in the tails. This is, again, due to the neglected higher-order terms, which are 
more relevant for higher values of the perturbation variables (i.e., in the tails of the pdf). Also 
in this case, the discrepancy is higher in the standard fiber because of the noise amplification 
that is possible with anomalous dispersion. 

10. Conclusions 

We have introduced a perturbative approach, based on the hydrodynamic formulation of the 
NLSE, to the nonlinear propagation of signal and noise in an optical fiber. In the small- 
dispersion approximation, the proposed method corresponds to the semiclassical approximation 
of the Madelung equations in quantum mechanics. However, we have seen that such approxi¬ 
mation is well-behaved only in the defocusing case (normal dispersion), for which the resulting 
first-order PDF system has a hyperbolic character. On the other hand, by removing this approx¬ 
imation, the method gets a more general validity and its accuracy is greatly improved. 

The developed formalism is rather general and, in principle, can be applied to input signals 
of any form. In this paper we have considered particular input data, namely the CW signal and 
the solitonic pulse, for which simple explicit equations have been derived. In the CW case, in 
particular, we have obtained explicit expressions for the PSD and pdf of the received signal 
that have been compared with numerical simulations (performed by employing the split-step 
Fourier method). The numerical results confirm the soundness of the proposed approach. 
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